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A fundamental prediction of relativistic cosmologies is that, owing to 

the expansion of space, observations of the distant cosmos should be 
time dilated and appear to run slower than events in the local universe. 
While observations of cosmological supernovae unambiguously show 
the expected redshift-dependent time dilation, this has not been the case 
for other distant sources. Here we present the identification of cosmic 
time dilation in a sample of 190 quasars monitored for over two decades 
in multiple wavebands by assessing various hypotheses through Bayesian 
analysis. This detection counters previous claims that observed quasar 
variability lacked the expected redshift-dependent time dilation. Hence, 
as well as dismissing the claim that the apparent lack of the redshift 
dependence of quasar variability represents a substantial challenge to 
the standard cosmological model, this analysis further indicates that the 
properties of quasars are consistent with them being truly cosmologically 


distant sources. 


A fundamental consequence of the relativistic picture of expanding 
space is cosmological time dilation, where events in the distant universe 
appear to run slowly compared with those in the local cosmos’ ’. While 
this time dilation has been unambiguously detected inthe light curves 
exhibited by cosmologically distant supernovae* *, the appearance 
of time dilation in other cosmic sources is less conclusive. For exam- 
ple, while examinations of the light curves of gamma-ray bursts have 
generally shown consistency with the expected cosmological signa- 
ture, uncertainties in the detailed emission mechanism and expected 
light-curve characteristics mean that this detection has not been defini- 
tive (for example, refs. 9-13). Furthermore, the role of the more recently 
discovered fast radio bursts“ as ‘standard clocks’ is similarly limited by 
knowledge of the physical processes driving the output”. 

Quasars have been known to be variable sources since their dis- 
covery inthe 1960s", with emission arising froma relativistic accretion 
disk orbiting a supermassive black hole”. However, it has been claimed 
that the variability displayed by quasars over a broad range of redshifts 
does not show the expected cosmological time dilation’*°. This has 
led to the suggestion that quasar variability is not intrinsic, but is due 
to microlensing owing to the presence of cosmologically distributed 
black holes”. Others have stated that this points to more fundamental 


issues with our cosmological ideas (for example, refs. 23-25), with 
even the suggestion that quasars are not cosmologically distant and 
that their observed redshifts are due to mechanisms other than the 
expansion of space. 

In 2012, a study of the variability characteristics of asample of 13 
quasars observed behind the Magellanic Clouds as part of the MACHO 
microlensing programme was suggestive of the expected (1+ z) time 
dilation dependence, where z is the quasar redshift”°. However, with 
the small sample and relatively short monitoring period, this result is 
inconclusive. Recently, a new sample of the variability properties of 
quasars was presented as part of the Dark Energy Survey and comprises 
190 quasars, covering the redshift range z = 0.2-4.0 (ref. 27). These 
are drawn from the sample of more than 100,000 spectroscopically 
identified quasars with absolute magnitude M < -22 in the 300 deg? 
Sloan Digital Sky Survey (SDSS) of Stripe 82. Published as part of the 
SDSS Data Release 7 quasar catalogue, the physical properties of these 
quasars are presented in ref. 28. This includes the bolometric luminos- 
ity, which was determined through spectral fitting and correction from 
composite spectral energy distributions”. These quasars were photo- 
metrically observed between 1998 and 2020, and so for more than two 
decades, through the combination of multiple epochs of exposures 


‘Sydney Institute for Astronomy, School of Physics, The University of Sydney, Sydney, New South Wales, Australia. 7Department of Statistics, 


The University of Auckland, Auckland, New Zealand. 


e-mail: geraint.lewis@sydney.edu.au 


Nature Astronomy 


Article 


https://doi.org/10.1038/s41550-023-02029-2 


Table 1 | The marginal likelihoods for the various hypotheses 
considered in this paper 


Hypothesis n log Z 2/Zmax 
Ho o -366.12 9.3x10° 
H, 1 -354.53 1 

H, Free -356.52 0.14 

H3 = -390.13 3.5x10° 
H, 2 -358.36 2.2x107 


As described in more detail in ‘Methodology’, these were calculated with the diffusive nested 
sampling approach DNest4°°. Here, Z represents the Bayesian evidence, with the maximum 
value across the various hypotheses being Z max- 


from SDSS, Panoramic Survey Telescope and Rapid Response System 
land the Dark Energy Survey, with additional follow-up monitoring 
with the Dark Energy Camera on the Blanco 4m. 

The total dataset consists of roughly 200 photometric observa- 
tions of each quasar in multiple bands, although the cadence of these 
observations is very uneven over the observing period. To account for 
this when calculating characteristic timescales of the quasar variabili- 
ties, ref. 27 adopted a Gaussian process regression (for example, ref. 30) 
to interpolate the photometric data and the associated uncertainties 
between the observations (details are given in Appendix A2 of their 
paper). Each quasar light curve in each band is represented as a damped 
random walk (DRW)*"”; this is found to be an accurate description of 
quasar variability with only a mild dependence on the physical prop- 
erties of the quasars”. Practically, this defines the covariance matrix 
of the Gaussian process that describes the variability. With this, the 
Gaussian process regression software, celerite™’, is used to determine 
the characteristic DRW timescale, as well as the 16th and 84th percen- 
tiles of the distribution. Armed with these bolometric luminosities 
and variability timescales drawn from the DRW analysis, the goal of 
this paper is to search for the signature of cosmological time dilation 
of these distant sources. 


Results 

In the following analysis, we consider the redshift dependence of 
time dilation to be of the form (1+ 2)", where zis the redshift of the 
source. Clearly, for the expected cosmological dependence, n =1, 
while n = O demonstrates no redshift dependence, representative 
of the claims from several previous studies of quasar samples (for 
example, ref. 20). To explore the various possibilities, several distinct 
hypotheses are explored. These are: Ho, nis fixed at zero, represent- 
ing no redshift dependence on the observed quasar timescales; H,, 
nis fixed at one, representing the expected redshift dependence of 
the cosmological time dilation; H,, n is treated as a free parameter; 
H, and H,, nis fixed at -1 and 2, respectively. The final two cases rep- 
resent extreme cases where additional influences, such as quasar 
evolution, may substantially influence the observed time variability 
of quasars. 

As outlined in ‘Methodology’, these differing hypotheses were 
compared through the calculation of the Bayesian evidence” for each 
situation under consideration, with the results of these calculations 
presented in Table 1; in assessing the ratio of Bayesian evidences, a 
factor of 10-100 is considered a strong favouring of one hypothesis 
over another, whereas greater than 100 is decisive”. One immediate 
conclusion is that the favoured hypothesis is H,, the case where n= 1, 
which represents the expected redshift dependence of the cosmologi- 
cal time dilation. This is significantly favoured over the alternative Ho, 
with an evidence ratio greater than 10°, which represents the situation 
where there is no redshift dependence on the observed timescales of 
cosmological variability. Furthermore, H, is significantly favoured over 
the two extreme cases, H, and Hy. 


P(n) 


Fig. 1| Posterior distribution of the redshift dependence for the time dilation. 
The posterior distribution ofn, where the redshift dependence of the observed time 
dilation is given by (1 + z)", for the Bayesian exploration of H,, where this index is 
treated asa free parameter in the analysis. From this distribution, n = 1.28+9:28, 
where the best fit value is taken as the median (50th percentile), while the 
uncertainties represent the 16th and 84th percentiles. This distribution was 


determined through an exploration of the posterior probability space with DNest4”°. 


The posterior distribution of the redshift dependence for the time 
dilation, n, specifically H}, where this is treated as a free parameter, is 
presented in Fig. 1. Reflecting the analysis presented previously, the 
value of nis clearly offset from zero, indicative of a redshift dependence 
of the observed timescale of variability over the quasar sample. This 
posterior distribution, which may be summarized using n = 1.28+?-8, 
is consistent with the expected cosmological dependence with n=1, 
and the presented analysis significantly favours the presence of cos- 
mological time dilation of the observed quasar variability. 


Methodology 

Probing the fundamental nature of our Universe often calls on standard 
rulers or candles to allow us to determine the influence of expansion 
on observable quantities. In hunting for cosmological time dilation, a 
standard clock with a measurable timescale is required. However, the 
challenge with objects such as quasars, and other cosmological sources 
such as gamma-ray bursts, is the complexity of the physical processes 
driving their variability. For quasars, where variability arises in the 
stochastic processes in the relativistic disk orbiting a supermassive 
black hole, the resultant luminosity fluctuations could potentially be 
dependent ona multitude of physical properties, including the mass 
of the central black hole, the degree of accretion and the wavelength 
of the observations. 

To address this, the sample of quasars under consideration here 
was split into anumber of subsamples of objects with similar intrinsic 
properties in terms of their bolometric luminosity and the rest wave- 
length of observations. The observations under consideration were 
taken in the (g, r, i) wavebands, and for the purpose of this study, the 
rest wavelength in each of the observed wavebands is defined to be 


_ 4,720A 
=~ 14+z 


_ 6,415A 


_ 7,835A 
> TO Te 


a A Tpz 


(1) 


where zis the redshift of the quasar under consideration and the numer- 
ical values are representative of the observed wavebands. 

The quasar subsamples are presented graphically in Fig. 2, which 
presents the rest wavelengths of each quasar in the (g, r, i) bands versus 
their bolometric luminosity. Each quasar has been colour-coded with 
its variability timescale, Tppw, assessed by fitting each observed light 
curve in each band with a DRW (see ref. 27 for more detail). Underly- 
ing the quasar sample are the regions of 12 subsamples under con- 
sideration in the salmon pink colour. These were chosen to have a 
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Fig. 2 | Quasar subsamples as a function of rest wavelength and bolometric 
luminosity. The entire quasar sample under consideration as a function of rest 
wavelength and bolometric luminosity, colour-coded with the DRW timescale, 
Tprw- The underlying rectangles in salmon pink represent the boundaries of 
the subsamples used in the analysis presented in this paper. Inset: the labelled 
numbers of the fields. 


Table 2 | The properties of the survey quasar subsamples 


Subsample AA (A) Alog,o N, Az Alog,oltorw 
(Lso/Lo) (days)] 

1 900 >1,900 46.73 47.2 37 1.6074.15 2.68 24.03 
2 1,900>2,900 46.7>47.2 if 0.8133.00 2.837411 

3 900 >1,900 46.2>46.7 74 1.5573.98  2.6124.23 
4 1,900>2,900 46.2>467 111 11123.01 2.49>4.42 
5 2,900>3,900 46.2>467 22 11121.70 3.03 >3.93 
6 900>1,900 45.7> 46.2 30 71483280 2707371 

7 1,900>2,900 45.7>46.2 101 0.68>2.80 2.55>3.92 
8 2,90073,900 45.7>46.2 58 0.607169 2.63>4.03 
9 3,90074,900 45.7>46.2 11 0.6071.00 2.66>3.84 
10 1,900>2,900 45.2>46.7 27 0.6371.45 2.3074.30 
11 2,90073,900 45.2>46.7 31 0.4771.45 2.27>4.30 
12 3,90074,900 45.2>46.7 20 047>70.98 23334.20 


The boundaries of the subsamples are given by AA, the rest wavelength, and Alog,o(Lga/Lo), 
the bolometric luminosity. The remaining columns give the number of quasar light curves in 
each subsample, N,,, as well as the redshift range of those quasars, Az, and timescale for the 
observed variability as given by treating this as a DRW, Alogjol Tppw(days)]. 


width in rest wavelength and bolometric luminosity of AA=1,000A 
and Alogi (Lgo/Lo) = 0.5. Note that the regions are continuous, 
with the top-left subsample spanning AA = 900 A> 1,900 A, and 
Alogio(Lgo/Lo) = 46.7 > 47.2; the details of the subsamples are given in 
Table 2. From Fig. 2, it is clear that these subsamples encompass the 
majority of quasars presented in this survey, and note that the combi- 
nation of the three wavebands means that each subsample of quasars 
contains a broader distribution of redshifts than ifthe wavebands were 
considered individually. Hence this combination provides a redshift 
lever arm that constrains the presence of cosmological time dilation 
in each subsample. We also note that a by-eye examination of Fig. 2 is 
suggestive of a gradient in the variability timescale over the sample. 
Given that in each subsample the quasars possess similar rest wave- 
lengths and bolometric luminosities, we make the assumption that they 


also possess the same characteristic intrinsic timescales, and hence 
any difference in timescale for a particular quasar subsample is due 
to the influence of cosmic time dilation and will show the appropriate 
dependence upon redshift. Of course, the physics of quasar variability 
is likely to depend on a number of factors, and so this assumption is 
considering that these quasars will exhibit similar variability properties 
inthe mean; we discuss this point again at the conclusion of this study. 
For each quasar subsample (labelled with k), we model the 
observed variability timescales (that is, logiol Zprw (days)]) as 


My = Cy + Nlogo +2) (2) 


where zis the redshift of the quasar and n is the power of the expected 
cosmological term, that is (1 + z)”. This model represents the variability 
with a different normalization term, C, for each wavelength-lumi- 
nosity bin, but demands the same cosmological dependence in terms 
of redshift. 

For the five distinct hypotheses considered, the normalization 
terms, C, were allowed to vary, and so for the cases where n is con- 
sidered to be a fixed value, this corresponds to an exploration of a 
12-dimensional posterior probability distribution. Physically, this 
situation reflects the situation where each wavelength-luminosity 
bin has a differing characteristic timescale, but a redshift dependence 
dependent upon the chosen value of n. For the remaining hypotheses, 
H,, where C, and n are treated as free parameters, this corresponds to 
a13-dimensional posterior probability distribution to be explored. 

To calculate the Bayesian evidence (also known as marginal 
likelihood) for each of these hypotheses, it is necessary to define a 
likelihood. It is important to note that the presented measurements 
and uncertainties of Tprw (in log,, space) are not symmetrical. Hence 
we represented the probability of each distribution of log,o[Tpaw 
(days)] as a skewed Gaussian, specifically scipy.stats.skewnorm in 
the numerical approach, which is written in Python (represented as 
SN). The 16th, 50th and 84th percentiles of this skewed probability 
distribution function (PDF) were fitted to the given values via a 
straight-forward optimization. This resulted in uncertainties of these 
percentiles of typically less than 2-3%. Hence we can define the log 
of the likelihood as 


logL=>) } SN - logPDFIM,,Am(Lmek) (3) 


k l=g,r,im=1.. -Nq 


where k sums over each of the subsample regions, Lover the observed 
wavebands and m over the number of quasars, N,, in the sample. Also, 
O,mare the parameters for the skewed Gaussian representing the prob- 
ability distribution for l0g,o[Tprw (days)] for each of the quasars in each 
of the waveband, and l, m e k implies that the quasar should only be 
considered if its rest-frame properties place it in subsample k. 

To explore the posterior probability distribution and calculate the 
Bayesian evidence (also knownas marginal likelihood), we employed 
diffusive nested sampling (DNest4)*°, a variant of the nested sampling 
technique”. This allows for correct posterior sampling and marginal 
likelihood estimation even in the case where the constrained prior 
distributions are difficult to sample from or explore with Markov chain 
Monte Carlo. For simplicity, we use uniform priors over the normaliza- 
tion parameters, C, between 1 and 5, and for H, where n is treated as 
a free parameter, a uniform distribution over nis adopted between -1 
and 3; the posterior distribution of n for this hypothesis is presented in 
Fig. 1. The normalizations, C, are well constrained and are reproduced 
for completeness in Fig. 3 for H2. 


Conclusions 

This paper presents the detection of the cosmological dependence of 
the time dilation in a recent sample of almost 200 quasars. These were 
monitored in multiple wavebands over atwo-decade period, allowing 
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46.7 > 47.2 46.2 > 46.7 
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45.7 > 46.2 C, | 45.2 > 45.7 


Fig. 3 | Posterior distributions for the normalization parameters. The 
posterior distributions for the normalization parameters, C,, for H, where nis 
treated as a free parameter. The bolometric luminosity range for each of the 
normalization parameters is given in the top left of each panel (Table 2). As 
with Fig. 1, these were the result of the sampling of the posterior probability 
distribution with DNest4. 


the determination of a characteristic timescale by treating the observed 
quasar variability as a DRW. 

Through an assessment of the Bayesian evidence, it was found that 
the hypothesis considering the expected (1 + z) cosmological depend- 
ence provides a significantly better description of the data than the 
case where there is no dependence on redshift. In considering the 
redshift dependence of quasar variability to be of the form (1+ 2)", 
where n is treated as a free parameter, the posterior distribution is 
found to be n = 1.28*9-28, again consistent with the expected cosmo- 
logical expansion of space. This detection of the cosmic expansion 
directly imprinted onto the variability of quasars further demonstrates 
that their observed properties are consistent with them being luminous 
and variable sources at cosmological distances, and counters previous 
claims that quasar variability is not intrinsic, but instead is due to 
external influences or non-standard physics. This has an immediate 
impact on various claims, such as the presence of a cosmologically 
dominant population of microlensing black holes (for example, 
refs. 18,22) or more esoteric ideas about the framework of the Uni- 
verse”, and is further evidence that we inhabit an expanding relativistic 
Universe. 

We do note that our result of n = 1.28*?8 could be consistent 
with an offset from the expected cosmological value of n = 1and could 
potentially indicate the presence of additional factors suchas an evolu- 
tion of quasars over cosmic time in addition to the time dilation due to 
cosmic expansion. Of course, we could imagine that quasar evolution 
over cosmic time could be responsible for the observed redshift 
dependence of the DRW timescale, but as we are considering similar 
quasars interms of the bolometric luminosity and observed rest wave- 
length, it would be a curious coincidence for this evolution to result in 
a(1+z) dependence to spoof cosmic expansion. Furthermore, if quasar 
evolution were solely responsible for the observed DRW properties, 
then the resulting lack of the expected cosmic time dilation would 
present a severe challenge to our cosmological model. However, it is 
important to note that there are some potential correlations of the 
DRW timescales with the inferred intrinsic properties of the quasars 
(for example, ref. 33), although these are not strong, and more exten- 
sive photometric datasets in terms of the number of quasars and the 


duration of their photometric light curves will be required to cleanly 
separate the influence of cosmic expansion from quasar evolution. 

In closing, we note that the lack of detection of the time dilation 
of quasar variability in previous studies is potentially due to the rela- 
tively small sample size in terms of the number of quasars under con- 
sideration”, or the cadence of data sampling and characterization 
of the quasar variability’’”°. Built on the observations of ref. 27, this 
present study has demonstrated that we are now in an epoch where 
we have observations of a sufficiently large number of quasars span- 
ning a broad range in redshifts, and observed over extended periods 
and with a cadence that overcomes their stochastic nature and results 
in an accurate characterization of their variability, yielding a robust 
determination of the imprint of cosmological expansion on their light 
curves. Furthermore, with upcoming programmes such as the Vera 
Rubin Observatory Legacy Survey of Space and Time, the number of 
quasars observed at high temporal cadence will rapidly increase and 
the measurement of cosmological time dilation, and potentially the 
influence of quasar evolution, will become readily observable (for 
example, ref. 32). 


Data availability 

The source data for this project are available at https://zenodo.org/ 
record/5842449#.YipOg-jMJPY, with the details of the available FITS 
tables presented in ref. 27. Note that a revised version of this catalogue 
was recently released due to an error in some rest-frame quantities. 
This revision does not impact any of the research presented in this 
paper. The software for this project is available at https://github.com/ 
eggplantbren/QuasarTimeDilation. 


Code availability 

This project made use of several publicly available software packages, 
especially DNest4* to undertake the exploration of the posterior prob- 
ability space and calculate the Bayesian evidence by integrating across 
this space. Further software packages employed include matplotlib”, 
numpy’’ and scipy“. Initial explorations of the posterior probability 
space were undertaken with emcee with corner plots prepared with 
corner”. The software employed as part of this project will be made 
available on reasonable request to the corresponding author. 
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